

In [1]:

    
from numpy.linalg import inv
import numpy as np
from scipy.linalg import eig

from sklearn.datasets import make_blobs
from sklearn.metrics import pairwise_distances

from diffmaps_util import k, diag

Diffusion Distance
A distance function between any two points based on the random walk on the graph [1].

Diffusion map
Low dimensional description of the data by the first few eigenvectors [1].



In [2]:

    
n = 3



In [3]:

    
X, y = make_blobs(n_samples=n, cluster_std=.1, centers=[[1,1]])
X









    Out[3]:





array([[ 1.0434619 ,  0.89997857],
       [ 1.0724103 ,  1.10392303],
       [ 1.1065059 ,  0.97652068]])

Define a pairwise similarity matrix between points...



In [4]:

    
L = k(X, .9)

and a diagonal normalization matrix $D_{i,i} = \sum_j L_{i,j}$



In [5]:

    
D = diag(L)

Matrix M
$M = D^{-1}L$



In [6]:

    
M = inv(D).dot(L)

The matrix M is adjoint to a symmetric matrix
$M_s = D^{1/2}MD^{-1/2}$

M and M_s share the same eigenvalues.
Since M_s is symmetric, it is diagonalizable and has a set of n real eigenvalues {$\lambda_{j=0}^{n-1}$} whose corresponding eigenvectors form an orthonormal basis of $\mathbf{R}^n$.
The left and right eigenvectors of M, denoted $\phi_j$ and $\psi_j$ are related to those of M_s.

$$ \phi_j = \mathbf{v}_j D^{1/2}, \psi_j = \mathbf{v}_j D^{-1/2} $$



In [7]:

    
Ms = diag(D, .5).dot(M).dot(diag(D,-.5))

Now we utilize the fact that by constrution M is a stochastic matrix



In [8]:

    
p0 = np.eye(n)

The stationary probability distribution $\Phi_0$



In [10]:

    
e = p0
for i in range(1000):
    e = e.dot(M)
print e









    



[[ 0.33309641  0.33309468  0.33380891]
 [ 0.33309641  0.33309468  0.33380891]
 [ 0.33309641  0.33309468  0.33380891]]



In [20]:

    
p1 = p0.dot(M)
p1









    Out[20]:





array([[ 0.33977989,  0.3241324 ,  0.33608772],
       [ 0.32504628,  0.34073789,  0.33421584],
       [ 0.33304245,  0.33025638,  0.33670117]])



In [15]:

    
w, v = eig(M)
w = w.real
print w
print v









    



[ 1.          0.01585329  0.00136565]
[[-0.57735027 -0.65564879 -0.48804677]
 [-0.57735027  0.7495651  -0.32879168]
 [-0.57735027 -0.090977    0.80852111]]



In [16]:

    
# sorting the eigenvalues and vectors
temp = {_:(w[_], v[:,_]) for _ in range(len(w))}
w = []
v = []
for _ in sorted(temp.items(), key=lambda x:x[1], reverse=True):
    w.append(_[1][0])
    v.append(_[1][1])
w = np.array(w)
v = np.array(v).T
print w
print v









    



[ 1.          0.01585329  0.00136565]
[[-0.57735027 -0.65564879 -0.48804677]
 [-0.57735027  0.7495651  -0.32879168]
 [-0.57735027 -0.090977    0.80852111]]



In [17]:

    
psi = v / v[:,0]
print psi









    



[[ 1.          1.13561702  0.8453218 ]
 [ 1.         -1.29828484  0.5694839 ]
 [ 1.          0.15757678 -1.40039963]]

Diffusion Map

$$ \Psi_t(x) = (\lambda_1^t\psi(x), \lambda_2^t\psi(x), ..., \lambda_k^t\psi(x)) $$



In [25]:

    
diffmap = (w.reshape(-1,1) * psi.T).T[:,1:]
diffmap









    Out[25]:





array([[ 0.01800327,  0.00115441],
       [-0.02058209,  0.00077771],
       [ 0.00249811, -0.00191245]])

Diffusion Distance

Defined by Euclidean distance in the diffusion map $$ D_t^2(x_0, x_1) = ||\Psi_t(x_0) - \Psi_t(x_1)||^2 $$



In [28]:

    
dt0 = pairwise_distances(diffmap)**2
dt0









    Out[28]:





array([[ 0.        ,  0.00148897,  0.00024982],
       [ 0.00148897,  0.        ,  0.00053993],
       [ 0.00024982,  0.00053993,  0.        ]])

Diffusion Distance [2]

Defined by probability distribution on time t. $$ D_t^2(x_0, x_1) = ||p(t, y|x_0) - p(t, y|x_1)||_w^2 \\ = \sum_y (p(t, y|x_0) - p(t, y|x_1))^2 w(y) $$



In [27]:

    
dt = []
for i in range(n):
    _ = []
    for j in range(n):
        _.append(sum((p1[i]-p1[j])**2 / v[:,0]**2))
    dt.append(_)
dt = np.array(dt)
dt









    Out[27]:





array([[ 0.        ,  0.00148898,  0.00024982],
       [ 0.00148898,  0.        ,  0.00053993],
       [ 0.00024982,  0.00053993,  0.        ]])



In [124]:

    
(dt0 - dt)









    Out[124]:





array([[  0.00000000e+00,  -6.69401358e-08,  -1.07584647e-08],
       [ -6.69401358e-08,   0.00000000e+00,  -2.21797167e-08],
       [ -1.07584647e-08,  -2.21797167e-08,   0.00000000e+00]])



In [147]:

    
print M
M.sum(axis=1)









    



[[ 0.34212144  0.32993451  0.32794404]
 [ 0.32540252  0.33742206  0.33717542]
 [ 0.3240756   0.33783864  0.33808576]]






    Out[147]:





array([ 1.,  1.,  1.])



In [153]:

    
w, v = eig(M)
w = w.real
print w
print v









    



[  1.00000000e+00   1.74591160e-02   1.70144248e-04]
[[ 0.57735027  0.81811916  0.05509647]
 [ 0.57735027 -0.35678015 -0.73203841]
 [ 0.57735027 -0.45098666  0.67903177]]



In [189]:

    
p0*w[0]*v[:,0]**2 + p0*w[1]*v[:,1]**2 + p0*w[2]*v[:,2]**2









    Out[189]:





array([[ 0.34501957,  0.        ,  0.        ],
       [ 0.        ,  0.33564692,  0.        ],
       [ 0.        ,  0.        ,  0.33696278]])